
// libMesh + MPI + PETSc includes
#include "moab_project_soln.h"

// iMesh + MOAB includes
#include "iMeshP.h"
#include "iMeshP_extensions.h"
#include "Coupler.hpp"
#include "MBiMesh.hpp"
#include "MBParallelConventions.h"

// libMesh includes
#include "vtk_io.h"
#include "gmsh_io.h"
#include "linear_implicit_system.h"
#include "parallel.h"

// libMesh and MOAB wrapper
#include "stacktrace.h"
#include "moab_libmesh_wrapper.h"

// C++ includes
#include <string>
#include <iostream>
#include <vector>
#include <algorithm>
#include <assert.h>
#include <sstream>


void MOABProjection::reduceMax(double &v){
  double buf;

  MPI_Barrier(PETSC_COMM_WORLD);
  MPI_Allreduce( &v, &buf, 1, MPI_DOUBLE, MPI_MAX, PETSC_COMM_WORLD);

  v=buf;
}


int MOABProjection::copy_from_to_moab_tag(moab::Core* moab, EquationSystems& es, Tag tagh, Range& entities,
            bool to_moab,
            std::map<int,int>& moab_libmesh_map)
{
	int result;

	// Get a reference to the LinearImplicitSystem we are solving
	LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>(PRIMARY_SYSTEM_NAME);

	PetscVector<PetscScalar>& src_soln = dynamic_cast<PetscVector<PetscScalar>&> (*system.solution) ;

	PetscInt rstart, rend, rank;
	result = VecGetOwnershipRange(src_soln.vec(),&rstart,&rend);CHKERRQ(result);

	MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

	std::vector<PetscScalar> ordered_vals(entities.size());

	//	PetscPrintf(PETSC_COMM_SELF, "\nProcessor %D -- Local size = %D; [ %D, %D ]", rank, lsize, rstart, rend) ;
	//	PetscPrintf(PETSC_COMM_SELF, "\nLocal size = %D; Range size = %D", lsize, entities.size()) ;

	if (to_moab)
	{
		system.update() ;

		std::vector<PetscInt> local_dofs(entities.size());
		std::vector<PetscScalar> sol_values;
		src_soln.localize(sol_values) ;

		for (unsigned int i=0; i < entities.size(); ++i)
		{
			std::map<int,int>::iterator it = moab_libmesh_map.find(entities[i]) ;
			if (it == moab_libmesh_map.end())
			{
				PetscPrintf(PETSC_COMM_SELF, "\n[%D] tag_set_data: Did not find key = %d", rank, entities[i]) ;
				return -1 ;
			}

			local_dofs[i] = it->second  ;
			ordered_vals[i] = sol_values[it->second];
		}

    // result = VecGetValues(src_local_soln.vec(), lsize, &local_dofs[0], &ordered_vals[0]) ;CHKERRQ(result);

		local_dofs.clear();

		// set the values from the source mesh solution
		result = moab->tag_set_data(tagh, entities, &ordered_vals[0]); PRINT_LAST_ERROR;
	}
	else
	{
		// get the values from the source mesh solution
		result = moab->tag_get_data(tagh, entities, &ordered_vals[0]); PRINT_LAST_ERROR;

//		int index = 0;
		for (unsigned int i=0; i < entities.size(); ++i)
		{
        std::map<int,int>::iterator it = moab_libmesh_map.find(entities[i]) ;
        if (it == moab_libmesh_map.end())
        {
          PetscPrintf(PETSC_COMM_SELF, "\n[%D] tag_get_data: Did not find key = %d", rank, entities[i]) ;
          return -1 ;
        }

        if (it->second >= rstart && it->second < rstart+rend )
        {
          // local_dofs[index++] = it->second  ;
          result = VecSetValue(src_soln.vec(), it->second, ordered_vals[i], INSERT_VALUES) ;CHKERRQ(result);
        }
		}

/*
		//assert(lsize == index) ;
    if (lsize != index)
    {
      PetscPrintf(PETSC_COMM_SELF, "\nProcessor %D -- Local size = %D; [ %D, %D ] but n(entities) = %D and index = %D", rank, lsize, rstart, rend, entities.size(), index) ;
  	  CHKERRQ(-1);
    }
*/

		// result = VecSetValues(src_soln.vec(), index, &local_dofs[0], &ordered_vals[0], INSERT_VALUES) ;CHKERRQ(result);

		// now update the current local solution
		src_soln.close() ;
		system.update() ;
	}

	ordered_vals.clear() ;

	return result ;

}


int MOABProjection::copy_from_to_moab_tag(moab::Core* moab, EquationSystems& es, std::string tag_name, Range& entities, bool to_moab, std::map<int,int>& moab_libmesh_map)
{
	int result;
	Tag tagh;

	result = moab->tag_get_handle(tag_name.c_str(), 1, MB_TYPE_DOUBLE, tagh); PRINT_LAST_ERROR;

	return copy_from_to_moab_tag(moab, es, tagh, entities, to_moab, moab_libmesh_map) ;
}


int MOABProjection::find_moab_libmesh_index_map(EquationSystems& es, moab::Core* moab, const Range& verts, std::map<int,int>& moab_libmesh_map)
{
	int result, rank;
	Tag global_id;

	this->projection_log->push("moab_libmesh_index_map()");

	MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

	moab_libmesh_map.clear();

	// Get a constant reference to the mesh object.
	const MeshBase& mesh = es.get_mesh();
	const PointLocatorBase& tree = mesh.point_locator();

	result = moab->tag_get_handle ("GLOBAL_ID", 1, moab::MB_TYPE_INTEGER, global_id, moab::MB_TAG_DENSE) ;PRINT_LAST_ERROR;

	std::vector<int> global_id_array(verts.size());
	result = moab->tag_get_data( global_id, verts, &global_id_array[0] ) ;

	std::vector<int>::iterator id_iter;

	std::vector<double> node_coords(verts.size() * 3) ;
	result = moab->get_coords(verts, &node_coords[0]);PRINT_LAST_ERROR;

	int search_method = 1;
	PetscOptionsGetInt("", "-index_search", &search_method,0);

	//PetscPrintf(PETSC_COMM_WORLD, "\nComputing projection between MOAB and libMesh\n");
	unsigned int index = 0 ;
	id_iter = global_id_array.begin();
	for (Range::iterator i = verts.begin(); i != verts.end(); ++i, ++id_iter, index+=3)
	{
		*id_iter = moab->id_from_handle( *i );

		//double node_coords[3]; // this will have the positions we are interested in
		//result = moab->get_coords(&(*i), 1, node_coords); PRINT_LAST_ERROR;

		int nid = -1;
		if (search_method == 1)
		{
			Point moab_point(node_coords[index], node_coords[index+1], node_coords[index+2]);

			const Elem* elem = tree(moab_point);

			internal_node_compare ndcompare(&node_coords[index]/*moab_point*/);
			for (unsigned int inode=0; inode < elem->n_nodes(); ++inode)
				if (ndcompare(elem->get_node(inode)))
					nid = elem->get_node(inode)->dof_number(0, 0, 0) ;
		}
		else
		{
			MeshBase::const_node_iterator       nd     = mesh.active_nodes_begin();
			const MeshBase::const_node_iterator end_nd = mesh.active_nodes_end();

			if (search_method == 2)
			{
				//MeshBase::const_node_iterator node_iter = std::find(mesh.active_nodes_begin(), mesh.active_nodes_end(), moab_point) ;
				//MeshBase::const_node_iterator node_iter = std::find(mesh.active_nodes_begin(), mesh.active_nodes_end(), [moab_point](Node const& n){return ((n-moab_point).size_sq() < tolerance);}) ;
				internal_node_compare ndcompare(&node_coords[index]/*moab_point*/);
				MeshBase::const_node_iterator node_iter = std::find_if(mesh.active_nodes_begin(), mesh.active_nodes_end(), ndcompare) ;

				if (node_iter != end_nd)
					nid = (*node_iter)->dof_number(0, 0, 0) ;
			}
			else
			{
				for ( ; nd != end_nd; ++nd)
				{
				  const Node* node = *nd;

				  if (fabs((*node)(0)-node_coords[index]) < 1e-5 && fabs((*node)(1)-node_coords[index+1]) < 1e-5
						  && fabs((*node)(2)-node_coords[index+2]) < 1e-5 )
					{
						nid = node->dof_number(0, 0, 0) ;
						break ;
					}
				}

			}
		}

		if (nid == -1) {
		  PetscPrintf(PETSC_COMM_SELF, "\n[%D] Couldn't find libMesh node matching MOAB node id = %D\t ", rank, *id_iter);
		}

		// use the MOAB global id as the key and the libMesh dof_number as value
		moab_libmesh_map.insert(std::pair<unsigned int, unsigned int>(*id_iter, nid));
		//libmesh_moab_map.insert(std::pair<unsigned int, unsigned int>(nid, *id_iter));

		//if (rank)
		//  PetscPrintf(PETSC_COMM_SELF, "\n[%D] Local Size = %D\tMOAB id = %D\t libMesh id = %D", rank, verts.size(), *id_iter, nid) ;
	}
	PetscBarrier(0);

	this->projection_log->pop("moab_libmesh_index_map()");

	return 0;
}


int MOABProjection::apply(std::string output_file)
{
	int result;
	int nprocs, rank;

	iMesh_Instance mesh;
	moab::Core* moab = new moab::Core();

	std::string readOpts, writeOpts;
	std::string ser_readOpts("");
	std::string par_readOpts(" moab:PARALLEL=READ_PART moab:PARTITION=PARALLEL_PARTITION moab:PARALLEL_RESOLVE_SHARED_ENTS moab:PARTITION_DISTRIBUTE moab:PARALLEL_GHOSTS=3.0.1 moab:CPUTIME ");
//		std::string par_readOpts(" moab:PARALLEL=READ_PART moab:PARTITION=PARALLEL_PARTITION moab:PARALLEL_RESOLVE_SHARED_ENTS moab:PARTITION_DISTRIBUTE  moab:CPUTIME DEBUG_IO=3 ");
	std::string par_writeOpts(" moab:PARALLEL=WRITE_PART moab:CPUTIME ");
	std::string ser_writeOpts("");
	std::string newReadOpt;

	const char *tmpFiles[2] = {"src.msh", "target.msh"};
	const char *inFiles[2] = {"temp1.h5m", "temp2.h5m"};
	const char *tag_name = "SOLUTION";
	const char *tag_options = "moab:TAG_STORAGE_TYPE=DENSE moab:TAG_DEFAULT_VALUE=0.0" ;

	PetscPrintf(PETSC_COMM_WORLD, "\n[Starting the MOAB mesh projection via MBCoupler]\n" );

	const unsigned int dim = source->get_mesh().mesh_dimension();

	MPI_Comm_size(PETSC_COMM_WORLD, &nprocs);
	MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
	std::ostringstream opt_str;
	opt_str << "DEBUG_IO= " << (verbose ? 3 : 0) << " " << (nprocs > 1 ? par_readOpts : ser_readOpts) ;
	readOpts = opt_str.str() ;
	opt_str.str("") ;
	opt_str << "DEBUG_IO= " << (verbose ? 3 : 0) << " " << (nprocs > 1 ? par_writeOpts : ser_writeOpts) ;
	writeOpts = opt_str.str() ;

	if (!use_new_method)
	{
		this->projection_log->push("write_libmesh_mesh()");

		// write out the source and target meshes so that we can read it back via MOAB iMesh interfaces
		// first convert to gmsh format and do mbconvert to h5m
		MeshBase& inFile = source->get_mesh();
		inFile.write(tmpFiles[0]);

		MeshBase& outFile = target->get_mesh();
		outFile.write(tmpFiles[1]);

		this->projection_log->pop("write_libmesh_mesh()");

		this->projection_log->push("convert_to_moab_mesh()");
		if (!rank)
		{
			std::vector<const char*> convert_args;

			// set commond options
			convert_args.push_back( "mbconvert" );
			//    convert_args.push_back( "-g" ) ;

			convert_args.push_back( tmpFiles[0] ) ;
			convert_args.push_back( inFiles[0] ) ;
			mbconvert(convert_args.size(), &convert_args[0]);

			convert_args[1] = tmpFiles[1];
			convert_args[2] = inFiles[1];
			mbconvert(convert_args.size(), &convert_args[0]);
		}
		this->projection_log->pop("convert_to_moab_mesh()");

		Parallel::barrier();
	}

	// read it in via the iMesh object
	mesh = reinterpret_cast<iMesh_Instance>( new MBiMesh(moab) );

	// array of file_sets and partition handles for source and target
	std::vector<EntityHandle> roots(2) ;
	std::vector<iMeshP_PartitionHandle> partn_handles(2);

	// array of parallel COMM instances for source and target
	MPI_Comm comms[2] = {PETSC_COMM_WORLD, PETSC_COMM_WORLD};
	std::vector<moab::ParallelComm *> pcs(2);

	MOAB_libMesh_Wrapper *moab_wrapper = NULL;
	if (use_new_method)
	{
		moab_wrapper = new MOAB_libMesh_Wrapper(moab, 2) ;

		std::vector<EntityHandle>& roots_wrapper = moab_wrapper->get_root_sets() ;
		std::vector<EntityHandle>& partition_handles = moab_wrapper->get_partition_handles() ;
		std::vector<moab::ParallelComm *>& parallel_communicators_wrapper = moab_wrapper->get_parallel_communicators() ;
		for (unsigned int index = 0; index < 2; index++) {
		  roots[index] = roots_wrapper[index];
		  partn_handles[index] = (iMeshP_PartitionHandle)partition_handles[index];
		  pcs[index] = parallel_communicators_wrapper[index];
		}
	}

	// load the meshes
	for (unsigned int index = 0; index < 2; index++) {

		if (use_new_method)
		{
			this->projection_log->push("convert_libmesh_to_moab");
			result = moab_wrapper->write_equation_systems_parallel(inFiles[index], (index == 0 ? *source : *target), index,
							true, verbose, verbose); PRINT_LAST_ERROR;
			this->projection_log->pop("convert_libmesh_to_moab");
		}
		else
		{
			// read in mesh(es)
			// Create root sets for each mesh using the iMesh API.  Then pass these
			// to the load_file functions to be populated.
			this->projection_log->push("read_moab_mesh");

			std::ostringstream extraOpt;
			extraOpt << readOpts ;

			// create the fileset handle
			result = moab->create_meshset( MESHSET_SET, roots[index] ); PRINT_LAST_ERROR;

			// create partitions and get corresponding ParallelComm objects;
			EntityHandle tmp_part;
			result = moab->create_meshset( MESHSET_SET, tmp_part ); PRINT_LAST_ERROR;
			pcs[index] = moab::ParallelComm::get_pcomm( moab, tmp_part, &comms[index] );
			partn_handles[index] = (iMeshP_PartitionHandle)tmp_part;

			if (nprocs > 1)
			{
				extraOpt << " moab:PCOMM=" << pcs[index]->get_id() ;
			}

			std::string newReadopts = extraOpt.str();

			// result = moab->create_meshset( moab::MESHSET_SET, file_sets[index]);PRINT_LAST_ERROR;
			PetscPrintf(PETSC_COMM_SELF, "\n[%D] [Reading mesh %s]: read_options : %s\n", rank, inFiles[index], newReadopts.c_str());

			iBase_EntitySetHandle file_set;
			file_set = (iBase_EntitySetHandle)roots[index];
			// Load the actual meshes
			iMesh_load(mesh, file_set, inFiles[index], newReadopts.c_str(), &result, strlen(inFiles[index]), newReadopts.length());
			if (result != 0) {
				PetscPrintf(PETSC_COMM_WORLD, "iMesh: Could not load the mesh file %s\n", inFiles[index]);
				PRINT_LAST_ERROR;
			}

			// perform a collective sync over all comm instances
			result = pcs[index]->collective_sync_partition();PRINT_LAST_ERROR;

			// assign the partition sets
			Range& part_set = pcs[index]->partition_sets();
			part_set.insert(roots[index]);

			// update the parallel comm instance about the local partitions set
			result = pcs[index]->set_partitioning(roots[index]);PRINT_LAST_ERROR;

			this->projection_log->pop("read_moab_mesh");
		}
	}

	PetscBarrier(0);

	// print details about the meshes
	PetscPrintf(PETSC_COMM_WORLD, "\nPrinting information about the entities\n");
	for (int index = 0; index < 2; index++) {
		moab::Range ents;
		// iterate over dimensions
		for (unsigned int d = 0; d <= dim; d++) {
			ents.clear();
			result = moab->get_entities_by_dimension((EntityHandle)roots[index], d, ents);PRINT_LAST_ERROR;
			PetscPrintf(PETSC_COMM_SELF, "\n%D\tFile %D\tFound %d %d-dimensional entities", rank, index, ents.size(), d);
			/*
			for (moab::Range::iterator it = ents.begin(); it != ents.end(); it++) {
			  moab::EntityHandle ent = *it;
			  PetscPrintf(PETSC_COMM_WORLD, "Found d=%d entity %d.\n", d, moab->id_from_handle(ent));
			}
			*/
		}
	}

	// add a tag named solution on both source and target meshes
	Tag tagh;
	bool created;
	double default_val = 0.0;
	result = moab->tag_get_handle (tag_name, 1, moab::MB_TYPE_DOUBLE, tagh, moab::MB_TAG_DENSE|moab::MB_TAG_CREAT, &default_val, &created) ;PRINT_LAST_ERROR;

	if (created) {
	  PetscPrintf(PETSC_COMM_WORLD, "\nCreating and adding a tag with name '%s' and of type double\n", tag_name);
	}
	else {
	  PetscPrintf(PETSC_COMM_WORLD, "\nTag with name '%s' already exists\n", tag_name);
	}

	this->projection_log->push("get_vertices()");
	// instantiate a coupler and load source elements into it
	moab::Range src_elems, src_verts, targ_elems, targ_verts, tmp_verts;

	result = pcs[0]->get_part_entities(src_elems, dim);PRINT_LAST_ERROR;
	result = moab->get_adjacencies(src_elems, 0, false, src_verts, Interface::UNION); PRINT_LAST_ERROR;
	//result = pcs[0]->filter_pstatus(src_verts, PSTATUS_NOT_OWNED, PSTATUS_NOT); PRINT_LAST_ERROR;

	result = pcs[1]->get_part_entities(targ_elems, dim);PRINT_LAST_ERROR;
	result = moab->get_adjacencies(targ_elems, 0, false, targ_verts, Interface::UNION); PRINT_LAST_ERROR;
	//result = pcs[1]->filter_pstatus(targ_verts, PSTATUS_NOT_OWNED, PSTATUS_NOT); PRINT_LAST_ERROR;

	PetscPrintf(PETSC_COMM_SELF, "\n[%D] Src Elements = %D and Target Elements = %D\n", rank, src_elems.size(), targ_elems.size());
	PetscPrintf(PETSC_COMM_SELF, "\n[%D] Src vertices = %D and Target vertices = %D\n", rank, src_verts.size(), targ_verts.size());
	this->projection_log->pop("get_vertices()");

	// find the index map between MOAB and libMesh representations
	// for both source and target meshes
	std::map<int,int> src_moab_libmesh_map;
	result = find_moab_libmesh_index_map(*source, moab, src_verts, src_moab_libmesh_map) ;

	std::map<int,int> targ_moab_libmesh_map;
	result = find_moab_libmesh_index_map(*target, moab, targ_verts, targ_moab_libmesh_map) ;

	// For each system, serialize the solution so we can project it onto
	// a potentially-very-different partitioning; we shall project only
	// the data from the first (implicit) system. Propagate the solution
	// from the ExplicitSystems as is since these will generally contain
	// only auxillary data.
	for (unsigned int i = 0; i < 1/*source->n_systems()*/; ++i)
	{
		System &old_sys = source->get_system(i);
		std::string current_sys_name = old_sys.name();

		PetscPrintf(PETSC_COMM_WORLD, "\nProjecting system '%s'\n", current_sys_name.c_str());
		PetscPrintf(PETSC_COMM_WORLD, " with variable %s\n", old_sys.variable_name(0).c_str());

		this->projection_log->push("tag_set_data()");

		// set the true solution from "src" to iMesh instance on the appropriate tag
		result = copy_from_to_moab_tag(moab, *source, tag_name, src_verts, true, src_moab_libmesh_map) ;CHKERRQ(result);

		// store empty data on the target
		target->get_system<LinearImplicitSystem>(current_sys_name).solution->zero();
		result = copy_from_to_moab_tag(moab, *target, tag_name, targ_verts, true, targ_moab_libmesh_map) ;CHKERRQ(result);

		// check if we zero out the solution and then copy from the tag data, if we get the exact same solution
		/*
		source->get_system<LinearImplicitSystem>(current_sys_name).solution->zero();
		ExodusII_IO (source->get_mesh()).write_equation_systems ("resource_old.exo", src);
		result = copy_from_to_moab_tag(moab, src, tag_name, src_verts, false, src_moab_libmesh_map) ;CHKERRQ(result);
		ExodusII_IO (source->get_mesh()).write_equation_systems ("resource_new.exo", src);
		*/

		this->projection_log->pop("tag_set_data()");

		// Calling sequence for mbcoupler_test interpolation rouinte
		std::vector<const char *> ssTagNameValues;
		double instant_time=0.0, pointloc_time=0.0, interp_time=0.0, gnorm_time=0.0, ssnorm_time=0.0;
		double toler = 5.e-10;
		std::string interpTag(tag_name);
		std::string gNormTag, ssNormTag;

		// test interpolation and global normalization and subset normalization
		result = test_interpolation(moab, dim, interpTag, gNormTag, ssNormTag,
								  ssTagNameValues, ssTagNameValues, roots, pcs,
								  instant_time, pointloc_time, interp_time,
								  gnorm_time, ssnorm_time, toler, verbose);
		PRINT_LAST_ERROR;

		this->projection_log->push("tag_get_data()");

		// get the values from the source mesh solution
		result = copy_from_to_moab_tag(moab, *target, tag_name, targ_verts, false, targ_moab_libmesh_map) ;CHKERRQ(result);

		this->projection_log->pop("tag_get_data()");

		reduceMax(instant_time);
		reduceMax(pointloc_time);
		reduceMax(interp_time);

		PetscPrintf(PETSC_COMM_WORLD, "\nMax time : %g %g %g [%g] (inst loc interp -- [total] %d procs )\n", instant_time,
			                pointloc_time, interp_time, instant_time+pointloc_time+interp_time, nprocs);

	}

  // output mesh
  if (!output_file.empty()) {
    this->projection_log->push("write_mesh()");
    Range partSets;
    // only save the target mesh
    partSets.insert(roots[1]);
    iMesh_save(mesh, (iBase_EntitySetHandle)roots[1], output_file.c_str(), writeOpts.c_str(), &result, output_file.length(), writeOpts.length()) ;PRINT_LAST_ERROR;

    //result = moab->write_file(output_file.c_str(), NULL, writeOpts.c_str(), partSets);PRINT_LAST_ERROR;
    PetscPrintf(PETSC_COMM_WORLD, "\nWrote target mesh and solution to %s\n", output_file.c_str());
    this->projection_log->pop("write_mesh()");
  }

  if (verbose)
  {
    this->projection_log->push("write_mesh()");
    iMesh_save(mesh, (iBase_EntitySetHandle)roots[0], "test1.h5m", writeOpts.c_str(),
                     &result, strlen("test1.h5m"), writeOpts.length()); PRINT_LAST_ERROR;
    this->projection_log->pop("write_mesh()");
    this->projection_log->push("write_mesh()");
    iMesh_save(mesh, (iBase_EntitySetHandle)roots[1], "test2.h5m", writeOpts.c_str(),
                     &result, strlen("test2.h5m"), writeOpts.length()); PRINT_LAST_ERROR;
    this->projection_log->pop("write_mesh()");
  }

  // !! ALL DONE !! Now clean up the memory.
  PetscPrintf(PETSC_COMM_SELF, "\n[%D] All done. Cleaning up the memory\n", rank);
  this->projection_log->push("cleanup()");

  if (use_new_method)
  {
    pcs.clear();
    delete moab_wrapper;
  }
  else
  {
    for (unsigned int i = 0; i < 2; i++)
      delete pcs[i];
  }

  // delete the mesh entirely
  result = moab->delete_mesh() ;PRINT_LAST_ERROR;

	roots.clear();
	partn_handles.clear();
  delete moab ;

  this->projection_log->pop("cleanup()");

  return 0;
}


ErrorCode report_iface_ents(Interface *mbImpl,
                              std::vector<ParallelComm *> &pcs,
                              const bool print_results)
{
  Range iface_ents[6];
  ErrorCode result = MB_SUCCESS;

  moab::Core* moab = dynamic_cast<moab::Core*>(mbImpl);

    // now figure out which vertices are shared
  for (unsigned int p = 0; p < pcs.size(); p++) {
    for (int i = 0; i < 4; i++) {
      result = pcs[p]->get_iface_entities(-1, i, iface_ents[i]);

      if (MB_SUCCESS != result) {
        std::cerr << "get_iface_entities returned error on proc "
                  << pcs[p]->proc_config().proc_rank() << "; message: " << std::endl;
        std::string last_error;
        mbImpl->get_last_error(last_error);
        if (last_error.empty()) std::cerr << "(none)" << std::endl;
        else std::cerr << last_error << std::endl;
      }
      if (0 != i) iface_ents[4].merge(iface_ents[i]);
    }
  }

    // report # iface entities
  result = mbImpl->get_adjacencies(iface_ents[4], 0, false, iface_ents[5],
                                   Interface::UNION);PRINT_LAST_ERROR;

  int rank;
  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

  if (print_results || iface_ents[0].size() != iface_ents[5].size()) {
    std::cerr << "Proc " << rank << " iface entities: " << std::endl;
    for (int i = 0; i < 4; i++)
      std::cerr << "    " << iface_ents[i].size() << " "
                << i << "d iface entities." << std::endl;
    std::cerr << "    (" << iface_ents[5].size()
              << " verts adj to other iface ents)" << std::endl;
  }

  return result;
}


ErrorCode MOABProjection::test_interpolation(Interface *mbImpl,
							 unsigned int dim,
                             std::string &interpTag,
                             std::string &gNormTag,
                             std::string &ssNormTag,
                             std::vector<const char *> &ssTagNames,
                             std::vector<const char *> &ssTagValues,
                             std::vector<EntityHandle> &roots,
                             std::vector<ParallelComm *> &pcs,
                             double &instant_time,
                             double &pointloc_time,
                             double &interp_time,
                             double &gnorm_time,
                             double &ssnorm_time,
                             double & toler,
                             bool verbose)
{
  ErrorCode result;
  int nprocs, rank;
  MPI_Comm_size(PETSC_COMM_WORLD, &nprocs);
  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

  moab::Core* moab = dynamic_cast<moab::Core*>(mbImpl);

  if (verbose)
  {
	  this->projection_log->push("print_iface_ents()");
	  result = report_iface_ents(mbImpl, pcs, true); PRINT_LAST_ERROR;
  	this->projection_log->pop("print_iface_ents()");
  }

  // source is 1st mesh, target is 2nd
  Range src_elems, targ_elems, targ_verts;

	this->projection_log->push("get_vertices()");
  result = pcs[0]->get_part_entities(src_elems, dim);PRINT_LAST_ERROR;
	this->projection_log->pop("get_vertices()");

  this->projection_log->push("coupler_init()");
  double start_time = MPI_Wtime();

  // instantiate a coupler, which also initializes the tree
  Coupler mbc(mbImpl, pcs[0], src_elems, 0);

  // initialize spectral elements, if they exist
  bool specSou=false, specTar = false;
  result =  mbc.initialize_spectral_elements(roots[0], (EntityHandle)roots[1], specSou, specTar);

  instant_time = MPI_Wtime();
  this->projection_log->pop("coupler_init()");

  this->projection_log->push("get_vertices()");
  // get points from the target mesh to interpolate
  // we have to treat differently the case when the target is a spectral mesh
  // in that case, the points of interest are the GL points, not the vertex nodes
  std::vector<double> vpos; // this will have the positions we are interested in
  int numPointsOfInterest = 0;
  if (!specTar)
  {// usual case
    Range tmp_verts;

    // first get all vertices adj to partition entities in target mesh
    result = pcs[1]->get_part_entities(targ_elems, dim);PRINT_LAST_ERROR;
    result = mbImpl->get_adjacencies(targ_elems, 0, false, targ_verts,
                     Interface::UNION);PRINT_LAST_ERROR;

    //result = pcs[1]->filter_pstatus(targ_verts, PSTATUS_NOT_OWNED, PSTATUS_NOT); PRINT_LAST_ERROR;

    // get position of these entities; these are the target points
    numPointsOfInterest = (int)targ_verts.size();
    vpos.resize(3*targ_verts.size());
    result = mbImpl->get_coords(targ_verts, &vpos[0]); PRINT_LAST_ERROR;
  }
  else
  {
    // in this case, the target mesh is spectral, we want values
    // interpolated on the GL positions; for each element, get the GL points, and construct CartVect!!!
    result = pcs[1]->get_part_entities(targ_elems, dim);PRINT_LAST_ERROR;

    result = mbc.get_gl_points_on_elements(targ_elems, vpos, numPointsOfInterest); PRINT_LAST_ERROR;
  }
  this->projection_log->pop("get_vertices()");

  PetscPrintf(PETSC_COMM_WORLD, "\nNumber of target elemnts = %D, and target vertices = %D", targ_elems.size(), targ_verts.size());

  PetscPrintf(PETSC_COMM_WORLD, "\nAbout to locate points with coordinates %D\n", vpos.size());

  this->projection_log->push("locate_points()");
    // locate those points in the source mesh
  result = mbc.locate_points(&vpos[0], numPointsOfInterest, 0, toler); PRINT_LAST_ERROR;
  this->projection_log->pop("locate_points()");

  pointloc_time = MPI_Wtime();

  PetscPrintf(PETSC_COMM_WORLD, "\nFinished locating points in %G time. Now interpolating\n", pointloc_time);

    // now interpolate tag onto target points
  std::vector<double> field(numPointsOfInterest);

  this->projection_log->push("interpolate_points()");

  //QuadraticHex
  //if(interpTag == "vertex_field" ||specSou ){
  //result = mbc.interpolate(Coupler::LINEAR_FE, interpTag, &field[0]);PRINT_LAST_ERROR;

  result = mbc.interpolate(Coupler::QUADRATIC_FE, interpTag, &field[0]);PRINT_LAST_ERROR;

  // }else if(interpTag == "element_field"){
  //  result = mbc.interpolate(Coupler::PLAIN_FE, interpTag, &field[0]);PRINT_LAST_ERROR;
  // }else {
  //  std::cout << "Using tag name to determine type of source field at the moment... Use either vertex_field or element_field\n";
  //  result = MB_FAILURE;
  //  PRINT_LAST_ERROR;
  // }

  this->projection_log->pop("interpolate_points()");

  interp_time = MPI_Wtime();

  iBase_EntitySetHandle file_set1 = (iBase_EntitySetHandle)roots[0];
  iBase_EntitySetHandle file_set2 = (iBase_EntitySetHandle)roots[1];

    // do global normalization if specified
  if (!gNormTag.empty()) {
    this->projection_log->push("normalize_mesh()");
    // Normalize the source mesh
    result = (ErrorCode)mbc.normalize_mesh(file_set1, gNormTag.c_str(), Coupler::VOLUME, 4);PRINT_LAST_ERROR;

    // Normalize the target mesh
    result = (ErrorCode)mbc.normalize_mesh(file_set2, gNormTag.c_str(), Coupler::VOLUME, 4);PRINT_LAST_ERROR;
    this->projection_log->pop("normalize_mesh()");
  }

  gnorm_time = MPI_Wtime();

    // do subset normalization if specified

  if (!ssNormTag.empty()) {
    this->projection_log->push("normalize_submesh()");
    result = (ErrorCode)mbc.normalize_subset(file_set1,
                               ssNormTag.c_str(),
                               &ssTagNames[0],
                               ssTagNames.size(),
                               &ssTagValues[0],
                               Coupler::VOLUME,
                               4);
    PRINT_LAST_ERROR;

    result = (ErrorCode)mbc.normalize_subset(file_set2,
                               ssNormTag.c_str(),
                               &ssTagNames[0],
                               ssTagNames.size(),
                               &ssTagValues[0],
                               Coupler::VOLUME,
                               4);
    PRINT_LAST_ERROR;
    this->projection_log->pop("normalize_submesh()");
  }

  ssnorm_time = MPI_Wtime();

  ssnorm_time -= gnorm_time;
  gnorm_time -= interp_time;
  interp_time -= pointloc_time;
  pointloc_time -= instant_time;
  instant_time -= start_time;

  this->projection_log->push("tag_set_data()");
    // set field values as tag on target vertices
  if (specSou)
  {
    // create a new tag for the values on the target
    Tag tag;
    std::string newtag = interpTag +"_TAR";
	PetscPrintf(PETSC_COMM_WORLD, "\nInterpolator: Creating spectral vertex source tag %s", newtag.c_str());
    result = mbImpl->tag_get_handle(newtag.c_str(), 1, MB_TYPE_DOUBLE, tag, MB_TAG_CREAT|MB_TAG_DENSE);
    result = mbImpl->tag_set_data(tag, targ_verts, &field[0]);
    PRINT_LAST_ERROR;
  }
  else
  {
    if (!specTar)
    {
      // use original tag
      Tag tag;

      PetscPrintf(PETSC_COMM_WORLD, "\nInterpolator: Getting linear source tag %s\n", interpTag.c_str());
      /*
      std::vector<Tag> tag_handles;
      result = moab->tag_get_tags(tag_handles);PRINT_LAST_ERROR;
      for (unsigned int i=0; i < tag_handles.size(); ++i)
      {
        std::string tagname;
        result = moab->tag_get_name(tag_handles[i], tagname);
        PetscPrintf(PETSC_COMM_WORLD, "\n\t Tag[%D] = %s", i, tagname);
      }
      */

      result = moab->tag_get_handle(interpTag.c_str(), 1, MB_TYPE_DOUBLE, tag); PRINT_LAST_ERROR;

      //PetscPrintf(PETSC_COMM_WORLD, "\nInterpolator: Target verts %D and field vals %D", targ_verts.size(), field.size());
      result = moab->tag_set_data(tag, targ_verts, &field[0]);PRINT_LAST_ERROR;
    }
    else
    {
      // we have the field values computed at each GL points, on each element
      // in the target mesh
      // we need to create a new tag, on elements, of size _ntot, to hold
      // all those values.
      // so it turns out we need ntot. maybe we can compute it from the
      // number of values computed, divided by number of elements
      int ntot = numPointsOfInterest/targ_elems.size();
      Tag tag;
      std::string newtag = interpTag +"_TAR";
	    PetscPrintf(PETSC_COMM_WORLD, "\nInterpolator: Getting spectral element target tag %s", interpTag.c_str());
      result = mbImpl->tag_get_handle(newtag.c_str(), ntot, MB_TYPE_DOUBLE, tag, MB_TAG_CREAT|MB_TAG_DENSE);
      result = mbImpl->tag_set_data(tag, targ_elems, &field[0]);
      PRINT_LAST_ERROR;

    }
  }
  this->projection_log->pop("tag_set_data()");

  // print out the first five interpolated field data for verification
  for (unsigned int i=0; i < (field.size() < 5 ? field.size() : 5); ++i)
    PetscPrintf(PETSC_COMM_SELF, "[%D]\t %s[%D] = %G\n", rank, interpTag.c_str(), i, field[i]);
  field.clear();

    // done
  return MB_SUCCESS;
}


